{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# LAVD detection and geometric detection\n\nNaive method to reproduce LAVD(Lagrangian-Averaged Vorticity deviation).\nIn the current example we didn't remove a mean vorticity.\n\nMethod are described here:\n\n    - Abernathey, Ryan, and George Haller. \"Transport by Lagrangian Vortices in the Eastern Pacific\",\n      Journal of Physical Oceanography 48, 3 (2018): 667-685, accessed Feb 16, 2021,\n      https://doi.org/10.1175/JPO-D-17-0102.1\n    - `Transport by Coherent Lagrangian Vortices`_,\n      R. Abernathey, Sinha A., Tarshish N., Liu T., Zhang C., Haller G., 2019,\n      Talk a t the Sources and Sinks of Ocean Mesoscale Eddy Energy CLIVAR Workshop\n\n    https://usclivar.org/sites/default/files/meetings/2019/presentations/Aberernathey_CLIVAR.pdf\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from datetime import datetime\n\nfrom matplotlib import pyplot as plt\nfrom numpy import arange, isnan, ma, meshgrid, zeros\n\nfrom py_eddy_tracker import start_logger\nfrom py_eddy_tracker.data import get_demo_path\nfrom py_eddy_tracker.dataset.grid import GridCollection, RegularGridDataset\nfrom py_eddy_tracker.gui import GUI_AXES\n\nstart_logger().setLevel(\"ERROR\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "class LAVDGrid(RegularGridDataset):\n    def init_speed_coef(self, uname=\"u\", vname=\"v\"):\n        \"\"\"Hack to be able to identify eddy with LAVD field\"\"\"\n        self._speed_ev = self.grid(\"lavd\")\n\n    @classmethod\n    def from_(cls, x, y, z):\n        z.mask += isnan(z.data)\n        datas = dict(lavd=z, lon=x, lat=y)\n        return cls.with_array(coordinates=(\"lon\", \"lat\"), datas=datas, centered=True)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def start_ax(title=\"\", dpi=90):\n    fig = plt.figure(figsize=(12, 5), dpi=dpi)\n    ax = fig.add_axes([0.05, 0.08, 0.9, 0.9], projection=GUI_AXES)\n    ax.set_xlim(-6, 36), ax.set_ylim(31, 45)\n    ax.set_title(title)\n    return fig, ax, ax.text(3, 32, \"\", fontsize=20)\n\n\ndef update_axes(ax, mappable=None):\n    ax.grid()\n    if mappable:\n        cb = plt.colorbar(\n            mappable,\n            cax=ax.figure.add_axes([0.05, 0.1, 0.9, 0.01]),\n            orientation=\"horizontal\",\n        )\n        cb.set_label(\"LAVD at initial position\")\n        return cb\n\n\nkw_lavd = dict(vmin=0, vmax=2e-5, cmap=\"viridis\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Load data cube of 3 month\nc = GridCollection.from_netcdf_cube(\n    get_demo_path(\"dt_med_allsat_phy_l4_2005T2.nc\"),\n    \"longitude\",\n    \"latitude\",\n    \"time\",\n    heigth=\"adt\",\n)\n\n# Add vorticity at each time step\nfor g in c:\n    u_y = g.compute_stencil(g.grid(\"u\"), vertical=True)\n    v_x = g.compute_stencil(g.grid(\"v\"))\n    g.vars[\"vort\"] = v_x - u_y"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Particles\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Time properties, for example with advection only 25 days\nnb_days, step_by_day = 25, 6\nnb_time = step_by_day * nb_days\nkw_p = dict(nb_step=1, time_step=86400 / step_by_day, u_name=\"u\", v_name=\"v\")\nt0 = 20236\nt0_grid = c[t0]\n# Geographic properties, we use a coarser resolution for time consuming reasons\nstep = 1 / 32.0\nx_g, y_g = arange(-6, 36, step), arange(30, 46, step)\nx0, y0 = meshgrid(x_g, y_g)\noriginal_shape = x0.shape\nx0, y0 = x0.reshape(-1), y0.reshape(-1)\n# Get all particles in defined area\nm = ~isnan(t0_grid.interp(\"vort\", x0, y0))\nx0, y0 = x0[m], y0[m]\nprint(f\"{x0.size} particles advected\")\n# Gridded mask\nm = m.reshape(original_shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## LAVD forward (dynamic field)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lavd = zeros(original_shape)\nlavd_ = lavd[m]\np = c.advect(x0.copy(), y0.copy(), t_init=t0, **kw_p)\nfor _ in range(nb_time):\n    t, x, y = p.__next__()\n    lavd_ += abs(c.interp(\"vort\", t / 86400.0, x, y))\nlavd[m] = lavd_ / nb_time\n# Put LAVD result in a standard py eddy tracker grid\nlavd_forward = LAVDGrid.from_(x_g, y_g, ma.array(lavd, mask=~m).T)\n# Display\nfig, ax, _ = start_ax(\"LAVD with a forward advection\")\nmappable = lavd_forward.display(ax, \"lavd\", **kw_lavd)\n_ = update_axes(ax, mappable)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## LAVD backward (dynamic field)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lavd = zeros(original_shape)\nlavd_ = lavd[m]\np = c.advect(x0.copy(), y0.copy(), t_init=t0, backward=True, **kw_p)\nfor i in range(nb_time):\n    t, x, y = p.__next__()\n    lavd_ += abs(c.interp(\"vort\", t / 86400.0, x, y))\nlavd[m] = lavd_ / nb_time\n# Put LAVD result in a standard py eddy tracker grid\nlavd_backward = LAVDGrid.from_(x_g, y_g, ma.array(lavd, mask=~m).T)\n# Display\nfig, ax, _ = start_ax(\"LAVD with a backward advection\")\nmappable = lavd_backward.display(ax, \"lavd\", **kw_lavd)\n_ = update_axes(ax, mappable)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## LAVD forward (static field)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lavd = zeros(original_shape)\nlavd_ = lavd[m]\np = t0_grid.advect(x0.copy(), y0.copy(), **kw_p)\nfor _ in range(nb_time):\n    x, y = p.__next__()\n    lavd_ += abs(t0_grid.interp(\"vort\", x, y))\nlavd[m] = lavd_ / nb_time\n# Put LAVD result in a standard py eddy tracker grid\nlavd_forward_static = LAVDGrid.from_(x_g, y_g, ma.array(lavd, mask=~m).T)\n# Display\nfig, ax, _ = start_ax(\"LAVD with a forward advection on a static velocity field\")\nmappable = lavd_forward_static.display(ax, \"lavd\", **kw_lavd)\n_ = update_axes(ax, mappable)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## LAVD backward (static field)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lavd = zeros(original_shape)\nlavd_ = lavd[m]\np = t0_grid.advect(x0.copy(), y0.copy(), backward=True, **kw_p)\nfor i in range(nb_time):\n    x, y = p.__next__()\n    lavd_ += abs(t0_grid.interp(\"vort\", x, y))\nlavd[m] = lavd_ / nb_time\n# Put LAVD result in a standard py eddy tracker grid\nlavd_backward_static = LAVDGrid.from_(x_g, y_g, ma.array(lavd, mask=~m).T)\n# Display\nfig, ax, _ = start_ax(\"LAVD with a backward advection on a static velocity field\")\nmappable = lavd_backward_static.display(ax, \"lavd\", **kw_lavd)\n_ = update_axes(ax, mappable)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Contour detection\nTo extract contour from LAVD grid, we will used method design for SSH, with some hacks and adapted options.\nIt will produce false amplitude and speed.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "kw_ident = dict(\n    force_speed_unit=\"m/s\",\n    force_height_unit=\"m\",\n    pixel_limit=(40, 200000),\n    date=datetime(2005, 5, 18),\n    uname=None,\n    vname=None,\n    grid_height=\"lavd\",\n    shape_error=70,\n    step=1e-6,\n)\nfig, ax, _ = start_ax(\"Detection of eddies with several method\")\nt0_grid.bessel_high_filter(\"adt\", 700)\na, c = t0_grid.eddy_identification(\n    \"adt\", \"u\", \"v\", kw_ident[\"date\"], step=0.002, shape_error=70\n)\nkw_ed = dict(ax=ax, intern=True, ref=-10)\na.filled(\n    facecolors=\"#FFEFCD\", label=\"Anticyclonic SSH detection {nb_obs} eddies\", **kw_ed\n)\nc.filled(facecolors=\"#DEDEDE\", label=\"Cyclonic SSH detection {nb_obs} eddies\", **kw_ed)\nkw_cont = dict(ax=ax, extern_only=True, ls=\"-\", ref=-10)\nforward, _ = lavd_forward.eddy_identification(**kw_ident)\nforward.display(label=\"LAVD forward {nb_obs} eddies\", color=\"g\", **kw_cont)\nbackward, _ = lavd_backward.eddy_identification(**kw_ident)\nbackward.display(label=\"LAVD backward {nb_obs} eddies\", color=\"r\", **kw_cont)\nforward, _ = lavd_forward_static.eddy_identification(**kw_ident)\nforward.display(label=\"LAVD forward static {nb_obs} eddies\", color=\"cyan\", **kw_cont)\nbackward, _ = lavd_backward_static.eddy_identification(**kw_ident)\nbackward.display(\n    label=\"LAVD backward static {nb_obs} eddies\", color=\"orange\", **kw_cont\n)\nax.legend()\nupdate_axes(ax)"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.10.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}